1 CONFIGURATIONS

# Set the working directory and tool paths on your local computer.
WD <- '/Users/Alec/motrpac/20210826_pass1c-umich'
# Set the gsutil path

knitr::opts_knit$set(root.dir=WD)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(cache = FALSE)
1.0.0.0.1 CSS Top Styling
write_css = FALSE
if(write_css){
  writeLines("td, th { padding : 3px } th { background-color:white; color:black; border:1px solid black; text-align:center } td {color:black; border:1px solid black; word-wrap:break-word; white-space:nowrap; overflow: hidden; text-overflow: ellipsis; max-width:300px; text-align:left}", con=file.path(normalizePath(WD), "style/style.css"))
}
1.0.0.0.2 Overview
1.0.0.0.3 Goals of Analysis
  • Examine data dimensions
  • Examine zero, negative, and missing Values
  • Impute missing values
  • Visualize sample-by-sample correlations
  • Identify outliers
  • Identify major causes of variance (drift, sex, control group, hour)
  • Visualize sample-by-feature heatmaps under different transformations
  • Examine sample median and sd distributions between each transformation
  • Transform data
  • Remove outlier samples
  • Save the processed data

2 Prepare Environment & Load Data

2.1 Setup the Environment

################################################################################
##### Resources and Dependencies ###############################################
################################################################################
# Whether to knit document and display data
knit_time = TRUE

# Load dependencies
pacs...man <- c("tidyverse","kableExtra","devtools","MotrpacBicQC","impute","glue",
                "rethinking")
for(pac in pacs...man){
  suppressWarnings(suppressPackageStartupMessages(library(pac, character.only = TRUE)))
}

#browseVignettes("MotrpacBicQC")
############################################################
##### Functions ############################################
############################################################

# Name functions
select <- dplyr::select
map <- purrr::map
desc <- dplyr::desc
arrange <- dplyr::arrange
melt <- reshape2::melt
mutate <- dplyr::mutate
glue <- glue::glue

# Global options
options(dplyr.print_max = 100)
options(scipen=10000)

# Colors
redblue100<-rgb(read.table(paste0(WD,'/colors/redblue100.txt'),sep='\t',row.names=1,header=T))

2.2 Log Variables (1)

ds <- 'pass1a'
site <- 'umich'
tech <- 'rppos'
TIS <- 'PLA'
tis <- 'pla'

2.3 Load Phenotype data

# Phenotype Data (1A)
#########################
pheno_file <- glue("{WD}/data/20201021_pass1a-06-pheno-viallabel_steep.txt")
pheno_df1a <- read.csv(pheno_file, header = T, sep = '\t')
pheno_df1a <- pheno_df1a %>%
  mutate(tissue = case_when(Specimen.Processing.sampletypedescription == 'Brown Adipose' ~ 'badi',
                            Specimen.Processing.sampletypedescription == 'EDTA Plasma' ~ 'pla',
                            Specimen.Processing.sampletypedescription == 'Gastrocnemius' ~ 'gas',
                            Specimen.Processing.sampletypedescription == 'Heart' ~ 'hrt',
                            Specimen.Processing.sampletypedescription == 'Kidney' ~ 'kid',
                            Specimen.Processing.sampletypedescription == 'Liver' ~ 'liv',
                            Specimen.Processing.sampletypedescription == 'Lung' ~ 'lun',
                            Specimen.Processing.sampletypedescription == 'White Adipose' ~ 'wadi',
                            Specimen.Processing.sampletypedescription == 'PaxGene RNA' ~ 'pax',
                            Specimen.Processing.sampletypedescription == 'Hippocampus' ~ 'hip',
                            Specimen.Processing.sampletypedescription == 'Cortex' ~ 'cor',
                            Specimen.Processing.sampletypedescription == 'Hypothalamus' ~ 'hyp',
                            Specimen.Processing.sampletypedescription == 'Vastus Lateralis' ~ 'vas',
                            Specimen.Processing.sampletypedescription == 'Tibia' ~ 'tib',
                            Specimen.Processing.sampletypedescription == 'Aorta' ~ 'aor',
                            Specimen.Processing.sampletypedescription == 'Small Intestine' ~ 'sma',
                            Specimen.Processing.sampletypedescription == 'Adrenals' ~ 'adr',
                            Specimen.Processing.sampletypedescription == 'Colon' ~ 'col',
                            Specimen.Processing.sampletypedescription == 'Spleen' ~ 'spl',
                            Specimen.Processing.sampletypedescription == 'Testes' ~ 'tes',
                            Specimen.Processing.sampletypedescription == 'Ovaries' ~ 'ova'))
pheno_df1a$viallabel <- as.character(pheno_df1a$viallabel)

2.4 Load Data (Kidney, Lung, Heart, Gastroc, Liver, Plasma, White Adipose)

Phenotypic Data Available: Kidney, Heart, Lung, Gastroc, Liver, Plasma, White Adipose Abundance Data not available for white adipose

2.4.1 Load the sample order files

# #created in 20210910_pass1a-umich-sample-annotation_steep.Rmd
# order_file <- glue("{WD}/data/20200910_pass1a1c-sample-order_steep.txt")
# sample.order<-read.delim(order_file,header=T, sep="\t")

# Load the prior pass1a data (takes a few minutes)
pass1a_nested_file <- glue("{WD}/../20200915_metabolomics-pass1a/data/20201010_pass1a-metabolomics-countdata-nested_steep.rds")
pass1a_df <- readRDS(pass1a_nested_file)

# Sample Order
pla_rppos_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rppos') %>%
  filter(TISSUE == 'plasma') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)
rm(pass1a_df)

2.4.2 Load the matrices as .RData files

# UM-rppos
load(file = glue("{WD}/data/UM-rppos/UM_rppos.0.RData"))
pla1a.0 <- pla_rppos_pass1a.0

2.5 Remove Internal Standards

# Remove internal standards
is <- colnames(pla1a.0)[grepl("istd", colnames(pla1a.0), ignore.case = TRUE)]
# Subset matrix
pla1a.0 <- pla1a.0[, !colnames(pla1a.0) %in% is]

2.5.1 Create Annotations

# Collect the distances from reference samples
tmp.iter1a <- pla_rppos_pass1a.0.order %>%
  mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
  mutate(drift = ifelse(sample_type == 'QC-DriftCorrection', 1, 0)) %>%
  arrange(sample_order)
tmp.iter1a$right_p <- 0
for(i in 1:nrow(tmp.iter1a)){
  # set t
  if(tmp.iter1a[i,'reference'] == 1){
    t = 0
  }else if(tmp.iter1a[i,'reference'] == 0){
    t = t + 1
  }
  tmp.iter1a[i,"right_p"] <- t
}
t=0
tmp.iter1a$right_p_d <- 0
for(i in 1:nrow(tmp.iter1a)){
  # set t
  if(tmp.iter1a[i,'drift'] == 1){
    t = 0
  }else if(tmp.iter1a[i,'drift'] == 0){
    t = t + 1
  }
  tmp.iter1a[i,"right_p_d"] <- t
}
t=0

tmp.iter1a <- tmp.iter1a %>%
  arrange(desc(sample_order))
tmp.iter1a$left_p <- 0
for(i in 1:nrow(tmp.iter1a)){
  # set t
  if(tmp.iter1a[i,'reference'] == 1){
    t = 0
  }else if(tmp.iter1a[i,'reference'] == 0){
    t = t + 1
  }
  tmp.iter1a[i,"left_p"] <- t
}
t=0
tmp.iter1a$left_p_d <- 0
for(i in 1:nrow(tmp.iter1a)){
  # set t
  if(tmp.iter1a[i,'drift'] == 1){
    t = 0
  }else if(tmp.iter1a[i,'drift'] == 0){
    t = t + 1
  }
  tmp.iter1a[i,"left_p_d"] <- t
}
t=0

tmp.iter1a <- tmp.iter1a %>%
  rowwise() %>%
  mutate(min_p = min(c(left_p,right_p) ,na.rm = TRUE)) %>%
  mutate(sum_p = sum(c(left_p,right_p) ,na.rm = TRUE)) %>%
  mutate(min_p_d = min(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
  mutate(sum_p_d = sum(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
  arrange(sample_order)
tmp.join1a <- tmp.iter1a %>%
  select(sample_id, sample_order, left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)

# Collect the sample order for test+ref samples
tmp.ref1a <- pla_rppos_pass1a.0.order %>%
  arrange(sample_order) %>% 
  mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
  mutate(drift = ifelse(grepl('Drift', sample_type), 1, 0)) %>%
  mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
  mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
                          grepl('0 hr', Key.anirandgroup) ~ 0,
                          grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
                          grepl('1 hr', Key.anirandgroup) ~ 1,
                          grepl('4 hr', Key.anirandgroup) ~ 4,
                          grepl('7 hr', Key.anirandgroup) ~ 7,
                          grepl('24 hr', Key.anirandgroup) ~ 24,
                          grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
  left_join(y = tmp.join1a) %>%
  select(sample_id,sample_order, Registration.sex, control, time, drift, reference,
         left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)
N2 <- tmp.ref1a[,1] %>% unlist()
miss1 <- N2[!N2 %in% row.names(pla1a.0)] # Verify all samples are in the pla1a.0 file
miss1
## named character(0)
# sample.order %>% filter(sample_id == "90750016606")
print('Vice Versa:')
## [1] "Vice Versa:"
miss2 <- row.names(pla1a.0)[!row.names(pla1a.0) %in% N2]
miss2
## character(0)
N2 <- N2[!N2 %in% c(miss1,miss2)] # TODO: investigate why samples are missing, continue for now
# Reorder pla1a.0 by run order
pla1a.0 <- pla1a.0[N2,]
all(N2 == row.names(pla1a.0)) # Must be true
## [1] TRUE
tmp.ref1a <- tmp.ref1a %>%
  filter(sample_id %in% N2) %>%
    select(sample_order, Registration.sex, control, time, drift, reference,
           left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>% as.matrix()

# Collect the sample order for test samples
options(digits = 14)
tmp.sample1a <- pla_rppos_pass1a.0.order %>%
  filter(substr(sample_id, 1, 1) == '9') %>%
  arrange(sample_order) %>% 
  mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
  mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
                          grepl('0 hr', Key.anirandgroup) ~ 0,
                          grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
                          grepl('1 hr', Key.anirandgroup) ~ 1,
                          grepl('4 hr', Key.anirandgroup) ~ 4,
                          grepl('7 hr', Key.anirandgroup) ~ 7,
                          grepl('24 hr', Key.anirandgroup) ~ 24,
                          grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
  left_join(y = tmp.join1a) %>%
  mutate(sample_id = str_replace_all(sample_id, pattern = '-', replacement = '')) %>%
  mutate(sample_id = as.numeric(sample_id)) %>%
  select(sample_id, sample_order, Registration.sex, control, time, 
         left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>%
  as.matrix()

N1 <- row.names(pla1a.0)[substr(row.names(pla1a.0), 1, 1) == '9']
tmp.sample1a <- tmp.sample1a[tmp.sample1a[,1] %in% N1,]
# Reorder pla1a.0.nr by run order
pla1a.0.nr <- pla1a.0[N1,]
tmp.sample1a[,1][!tmp.sample1a[,1] %in% row.names(pla1a.0.nr)] # Verify all samples are in the pla1a.0 file
## numeric(0)
all(as.character(tmp.sample1a[,1]) == row.names(pla1a.0.nr))
## [1] TRUE
# If out of order, command below will ensure pla1a.0.nr in run order
#pla1a.0.nr <- pla1a.0.nr[as.character(tmp.sample1a[,1]),]

3 Dimensions, Zero/Neg/Missing Values, & Log2

3.1 Dimensions (with reference samples)

NR <- dim(pla1a.0)[1]
P <- dim(pla1a.0)[2]
dim(pla1a.0)
## [1] 101 170

3.2 Dimensions (without reference samples)

N <- dim(pla1a.0.nr)[1]
dim(pla1a.0.nr)
## [1]  77 170

3.3 Negative or Zero Values

confirmed: no negative or zero values

min(pla1a.0,na.rm=T)
## [1] 59.47415952

3.4 Missing Features (with references)

Blank reference samples at the beginning and end

pla1a.0.f.c0<-apply(pla1a.0,1,function(x) sum(is.na(x))) 
plot(pla1a.0.f.c0, ylim = c(0,P))

3.5 Blank Samples (without references)

No blank test samples

pla1a.0.nr.f.c0<-apply(pla1a.0.nr,1,function(x) sum(is.na(x))) 
plot(pla1a.0.nr.f.c0, ylim = c(0,P))

3.6 Examine Distribution of missing values (with reference samples)

pla1a.0.f.c0<-apply(pla1a.0,2,function(x) sum(is.na(x))) 
plot(pla1a.0.f.c0, ylim = c(0,NR))

3.7 Examine Distribution of missing values (test samples only)

pla1a.0.nr.f.c0<-apply(pla1a.0.nr,2,function(x) sum(is.na(x))) 
plot(pla1a.0.nr.f.c0, ylim = c(0,N)); abline(h = 10)

3.8 Remove high-missing features

Remove high missing features

rm_n <- sum(pla1a.0.nr.f.c0>=10)
pla1a.0 <- pla1a.0[,pla1a.0.nr.f.c0<10]
pla1a.0.nr <- pla1a.0.nr[,pla1a.0.nr.f.c0<10]
dim(pla1a.0)
## [1] 101 170
dim(pla1a.0.nr)
## [1]  77 170

3.9 Take the log2

pla1a.0.1 <- log(pla1a.0, 2)
pla1a.0.nr1 <- log(pla1a.0.nr, 2)

3.10 Total missing values (includes reference samples)

sum(is.na(pla1a.0.1))
## [1] 0

3.11 Total missing values (test samples)

sum(is.na(pla1a.0.nr1))
## [1] 0
feature_impute <- apply(is.na(pla1a.0.nr1),2,sum)[apply(is.na(pla1a.0.nr1),2,sum) > 0]

3.12 Confirm New Distribution of missing features (test samples only)

pla1a.0.nr1.f.c0<-apply(pla1a.0.nr1,2,function(x) sum(is.na(x))) 
plot(pla1a.0.nr1.f.c0, ylim = c(0,NR))

3.13 Log Variables (2)

# NR
# N
# P
neg_vals <- 0
zero_vals <- 0
feature_na_filter <- 10
P1 <- dim(pla1a.0.nr)[2]
NR1 <- dim(pla1a.0)[1]
N1 <- dim(pla1a.0.nr)[1]
na_vals_impute <- sum(is.na(pla1a.0.nr1))
knn_k <- 10
feature_impute
## named integer(0)

4 Imputation

4.1 Imputation (test samples)

if(na_vals_impute > 0){
  glue("Features & Values to impute:")
  feature_impute
  pla1a.0.nr2<-impute.knn(pla1a.0.nr1,k=10)$data
  #view the features before and after imputation
  par(mfrow=c(2,1),bg="black")
  image(as.matrix(pla1a.0.nr1[,pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
  image(as.matrix(pla1a.0.nr2[, pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
  par(mfrow=c(1,1) ,bg="white")
  glue("Verified no missing values: {sum(is.na(pla1a.0.nr2))}")  #verified 0
}else{
  print('No missing values to impute')
  pla1a.0.nr2 <- pla1a.0.nr1
}
## [1] "No missing values to impute"

4.2 Imputation (+ reference samples)

# feature_impute2 <- apply(is.na(pla1a.0.1),2,sum)[apply(is.na(pla1a.0.1),2,sum) > 0]
# glue("Features & Values to impute:")
# feature_impute2
# pla1a.0.2<-impute.knn(pla1a.0.1,k=10)$data
# #view the features before and after imputation
# par(mfrow=c(2,1),bg="black")
# image(as.matrix(pla1a.0.1[,pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
# image(as.matrix(pla1a.0.2[,pla1a.0.nr1.f.c0>0]),col=redblue100,axes=F)
# par(mfrow=c(1,1) ,bg="white")
# glue("Verified no missing values: {sum(is.na(pla1a.0.2))}")  #verified 0

5 NxN heatmaps

5.1 Plotting NxN heatmaps (+ reference samples)

a <- 0.92
b <- 1
x <- tmp.ref1a[,1]
sidebar  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar <- cbind(sidebar,sidebar,sidebar)
cor.tmp<-cor(t(pla1a.0.1),method="spearman",use="pairwise.complete.obs") #includes ref samples
dim(cor.tmp)
## [1] 101 101
image(cbind(cor.tmp,sidebar),
            col=redblue100,axes=FALSE,zlim=c(0.92,1), main=glue("{TIS}2-{NR1}, run-order, z=(0.92,1)"), asp = 1)

if(knit_time){
  pla_rppos_pass1a.0.order %>%
    arrange(sample_order) %>%
    select(sample_id, sample_type, sample_order) %>%
    knitr::kable(format = "html") %>%
    scroll_box(width = "100%", height = "400px")
}
sample_id sample_type sample_order
CS00000MP-1 QC-DriftCorrection 3
R00CHRPL2-8 QC-InternalStandard 4
CSMR80015-1 QC-Reference 5
CSMR80014-1 QC-Reference 6
CS00000MP-2 QC-DriftCorrection 7
CS000QCMP-1 QC-Pooled 8
90109013104 Sample 9
90121013104 Sample 10
90156013104 Sample 11
90031013104 Sample 12
90127013104 Sample 13
90133013104 Sample 14
90044013104 Sample 15
90025013104 Sample 16
90138013104 Sample 17
CS00000MP-3 QC-DriftCorrection 18
90159013104 Sample 19
90148013104 Sample 20
90114013104 Sample 21
90041013104 Sample 22
90005013104 Sample 23
90017013104 Sample 24
90050013104 Sample 25
90037013104 Sample 26
90142013104 Sample 27
CS00000MP-4 QC-DriftCorrection 28
CS000QCMP-2 QC-Pooled 29
90010013104 Sample 30
90027013104 Sample 31
90009013104 Sample 32
90123013104 Sample 33
90034013104 Sample 34
90124013104 Sample 35
90008013104 Sample 37
90117013104 Sample 38
CS00000MP-5 QC-DriftCorrection 39
90007013104 Sample 40
90112013104 Sample 41
90028013104 Sample 42
90011013104 Sample 43
90026013104 Sample 44
90141013104 Sample 45
90126013104 Sample 46
90115013104 Sample 47
90134013104 Sample 48
CS00000MP-6 QC-DriftCorrection 49
CS000QCMP-3 QC-Pooled 50
90042013104 Sample 51
90119013104 Sample 52
90023013104 Sample 53
90129013104 Sample 54
90029013104 Sample 55
90020013104 Sample 56
90046013104 Sample 57
90122013104 Sample 58
90039013104 Sample 59
CS00000MP-7 QC-DriftCorrection 60
90157013104 Sample 61
90048013104 Sample 62
90032013104 Sample 63
90018013104 Sample 64
90160013104 Sample 65
90155013104 Sample 66
90145013104 Sample 67
90140013104 Sample 68
90033013104 Sample 69
CS00000MP-8 QC-DriftCorrection 70
CS000QCMP-4 QC-Pooled 71
90110013104 Sample 72
90143013104 Sample 73
90147013104 Sample 74
90038013104 Sample 75
90139013104 Sample 76
90150013104 Sample 77
90053013104 Sample 78
90014013104 Sample 79
90052013104 Sample 80
CS00000MP-9 QC-DriftCorrection 81
90146013104 Sample 82
90040013104 Sample 83
90118013104 Sample 84
90120013104 Sample 85
90128013104 Sample 86
90043013104 Sample 87
90136013104 Sample 88
90013013104 Sample 89
90135013104 Sample 90
CS00000MP-10 QC-DriftCorrection 91
CS000QCMP-5 QC-Pooled 92
90015013104 Sample 93
90047013104 Sample 94
90152013104 Sample 95
90144013104 Sample 96
90012013104 Sample 97
90045013104 Sample 98
CS000QCMP-6 QC-Pooled 99
CS00000MP-11 QC-DriftCorrection 100
CSMR80015-2 QC-Reference 101
CSMR80014-2 QC-Reference 102
R00CHRPL1-3 QC-InternalStandard 103
CS00000MP-12 QC-DriftCorrection 106

5.2 Plotting NxN heatmaps (test samples)

a <- 0.92
b <- 1
x <- tmp.sample1a[,2]
sidebar  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar <- cbind(sidebar,sidebar,sidebar)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
glue("Verified {N} test samples: {dim(cor.tmp)}") #verified, =78
## Verified 77 test samples: 77
## Verified 77 test samples: 77
image(
  cbind(cor.tmp[order(tmp.sample1a[,2]),],sidebar),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Run-Order, z=(0.92,1)"), asp = 1)

plot(apply(cor.tmp[order(tmp.sample1a[,2]),order(tmp.sample1a[,2])],1,median))

5.3 Examine Outlier Samples (test samples)

cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
dim(cor.tmp) #verified, =78
## [1] 77 77
plot(apply(cor.tmp,1,median)); abline(h=0.96, col = 'blue')

# Determine which samples are outliers
o.n <- c(1:N1)[apply(cor.tmp,1,median)<0.96]
o.s <- colnames(cor.tmp)[apply(cor.tmp,1,median)<0.96]
glue("Outlier Samples: {length(o.s)}")
## Outlier Samples: 1
o.s
## [1] "90122013104"
if(knit_time){
  pla_rppos_pass1a.0.order %>% 
    filter(sample_id %in% o.s) %>%
    knitr::kable(format = "html") %>%
    scroll_box(width = "100%", height = "100%")
}
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90122013104 Sample 58 Exercise - 4 hr 1 
o.f <- 0.96

5.4 Visualize the Sex-Time-Control NxN Heatmap

The one sample that could be classified as an outlier seems to have a similar expression signature of is male/female contemporaries

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,3]
sex.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- tmp.sample1a[,4]
control.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- as.factor(tmp.sample1a[,5]) %>% as.numeric()
time.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(sex.type,sex.type,sex.type, 
               control.type,control.type,control.type, 
               time.type,time.type,time.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),
          order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5])],
  sidebar[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N}, Sex-Control-Time, z=(0.92,1)"), asp = 1)

plot(apply(cor.tmp[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]), 
                         order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5])],
           1,median))

5.5 NxN Heatmap: Test and Reference Samples by Run Order

# Add color sidebar
a <- 0.92
b <- 1
x <- is.na(tmp.ref1a[,2]) %>% as.integer()
sample.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(sample.type,sample.type,sample.type,sample.type,sample.type)
cor.tmp<-cor(t(pla1a.0.1),method="spearman",use="pairwise.complete.obs") #includes ref samples
image(
  cbind(
  cor.tmp, sidebar),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2, Run-Order (+Ref-type), z=(0.92,1)"), asp = 1)

5.6 Log Variables (3)

run_var <- 0
outlier_sample_n <- length(o.s)
outlier_samples <- o.s
outlier_filter <- o.f

6 Positioning of Experimental Samples in Relation to Reference/Drift Samples

6.1 Position Adjacent to Reference Sample (Left)

# Add color sidebar
a <- 0.92
b <- 1
x <- tmp.sample1a[,6]
left.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(left.type, left.type, left.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,6]),
          order(tmp.sample1a[,6])],
  sidebar[order(tmp.sample1a[,6]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Left Pos, z=(0.92,1)"), asp = 1)

6.2 Position Adjacent to Reference Sample (Right)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,7]
right.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(right.type, right.type, right.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,7]),
          order(tmp.sample1a[,7])],
  sidebar[order(tmp.sample1a[,7]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Right Pos, z=(0.92,1)"), asp = 1)

6.3 Position Adjacent to Reference Sample (Both-Min)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,8]
min.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(min.type, min.type, min.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,8]),
          order(tmp.sample1a[,8])],
  sidebar[order(tmp.sample1a[,8]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Min Pos, z=(0.92,1)"), asp = 1)

# Boxplot
plot(apply(cor.tmp[order(tmp.sample1a[,8]),order(tmp.sample1a[,8])],1,median))

x <- apply(cor.tmp[order(tmp.sample1a[,8]),order(tmp.sample1a[,8])],1,median)

df <- data.frame(sample_id = names(x), cor_median = x, min_pos = sort(as.integer(tmp.sample1a[,8]))) %>%
  left_join(y = pla_rppos_pass1a.0.order) %>%
  mutate(min_pos = factor(min_pos))
df %>% ggplot(aes(x = min_pos, y = cor_median)) +
  geom_boxplot(alpha = 0.5) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.8, color = 'forestgreen')

6.4 Position Adjacent to Reference Sample (Both-Sum AKA Block Length)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,9]
sum.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(sum.type, sum.type, sum.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,9]),
          order(tmp.sample1a[,9])],
  sidebar[order(tmp.sample1a[,9]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Block Length, z=(0.92,1)"), asp = 1)

6.5 Position Adjacent to Drift Sample (Left)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,10]
left.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(left.type, left.type, left.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,10]),
          order(tmp.sample1a[,10])],
  sidebar[order(tmp.sample1a[,10]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Left Pos (Drift), z=(0.92,1)"), asp = 1)

plot(apply(cor.tmp[order(tmp.sample1a[,10]),order(tmp.sample1a[,10])],1,median))

plot(apply(cor.tmp[order(tmp.sample1a[,10]),order(tmp.sample1a[,10])],1,median), ylim = c(0.925,1))

# Boxplot
x <- apply(cor.tmp[order(tmp.sample1a[,10]),order(tmp.sample1a[,10])],1,median)
df <- data.frame(sample_id = names(x), cor_median = x, left_pos = sort(as.integer(tmp.sample1a[,10]))) %>%
  left_join(y = pla_rppos_pass1a.0.order) %>%
  mutate(left_pos = factor(left_pos))
df %>% ggplot(aes(x = left_pos, y = cor_median)) +
  geom_boxplot(alpha = 0.5) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.8, color = 'forestgreen')

6.6 Position Adjacent to Drift Sample (Right)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,11]
right.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(right.type, right.type, right.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,11]),
          order(tmp.sample1a[,11])],
  sidebar[order(tmp.sample1a[,11]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Right Pos (Drift), z=(0.92,1)"), asp = 1)

6.7 Position Adjacent to Drift Sample (Both-Min)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,12]
min.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(min.type, min.type, min.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,12]),
          order(tmp.sample1a[,12])],
  sidebar[order(tmp.sample1a[,12]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Min Pos (Drift), z=(0.92,1)"), asp = 1)

6.8 Position Adjacent to Drift Sample (Both-Sum AKA Block Length)

# Add color sidebar for sex
a <- 0.92
b <- 1
x <- tmp.sample1a[,13]
sum.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(sum.type, sum.type, sum.type)
cor.tmp<-cor(t(pla1a.0.nr1),method="spearman",use="pairwise.complete.obs")
image(
  cbind(
  cor.tmp[order(tmp.sample1a[,13]),
          order(tmp.sample1a[,13])],
  sidebar[order(tmp.sample1a[,13]),]),
  col=redblue100,axes=FALSE, zlim=c(0.92,1), main=glue("{TIS}2-{N1}, Block Length (Drift), z=(0.92,1)"), asp = 1)

plot(apply(cor.tmp[order(tmp.sample1a[,13]),order(tmp.sample1a[,13])],1,median))

# Boxplot
x <- apply(cor.tmp[order(tmp.sample1a[,13]),order(tmp.sample1a[,13])],1,median)
df <- data.frame(sample_id = names(x), cor_median = x, drift_block_length = sort(as.integer(tmp.sample1a[,13]))) %>%
  left_join(y = pla_rppos_pass1a.0.order) %>%
  mutate(drift_block_length = factor(drift_block_length))
df %>% ggplot(aes(x = drift_block_length, y = cor_median)) +
  geom_boxplot(alpha = 0.5) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.8, color = 'forestgreen')

7 N-P Heatmaps, Median and SD Distributions, & Transformations

7.0.1 N-P Heatmaps: Log2, imputed, sample order,

Note: Samples as columns Sidebar: p-values, t-values, and sig (binary) comparing normal test samples and outlier test samples

n.s <- row.names(pla1a.0.nr2)[!row.names(pla1a.0.nr2) %in% o.s]
# dim(pla1a.0.nr2[o.s,])
# dim(pla1a.0.nr2[n.s,])
hmo <- heatmap(pla1a.0.nr2)$colInd

if(length(o.s) > 1){
  i <- 1
df_t <- data.frame()
for(i in 1:ncol(pla1a.0.nr2)){
  feat <- colnames(pla1a.0.nr2)[i]
  tobj <- t.test(y = pla1a.0.nr2[n.s,i], x = pla1a.0.nr2[o.s,i])
  df_t <- rbind(df_t, data.frame(feat, t_val = round(tobj$statistic,digits = 2),p_val = round(tobj$p.value, digits = 16)))
}
df_t <- df_t %>% arrange(t_val) %>%
  mutate(sig = ifelse(p_val <= 0.05, 1, 0))
row.names(df_t) <- df_t$feat

# Add color sidebar for pval
a <- range(pla1a.0.nr2)[1]
b <- range(pla1a.0.nr2)[2]
x <- df_t[colnames(pla1a.0.nr2)[hmo], "p_val"]
p.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- df_t[colnames(pla1a.0.nr2)[hmo], "t_val"]
t.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- df_t[colnames(pla1a.0.nr2)[hmo], "sig"]
sig.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-rbind(p.type, p.type,
               t.type, t.type,
               sig.type, sig.type)

image(rbind(pla1a.0.nr2[order(tmp.sample1a[,2]),hmo],sidebar)
      ,col=redblue100,axes=F,main=glue("{TIS}2, f{P1}-HC, Run-Order"), asp = 1)
}else{
  image(pla1a.0.nr2[order(tmp.sample1a[,2]),hmo],col=redblue100,axes=F,main=glue("{TIS}2, f{P1}-HC, Run-Order"), asp = 1)
}

plot(apply(pla1a.0.nr2[order(tmp.sample1a[,2]),hmo],1,median))

7.0.2 N-P Heatmap: Log2, imputed, Sex-Control-Time order (2)

Note: pla1a.0.2 and pla1a.0.nr2 represent knn-imputed and log2 Note: Samples as columns

a <- range(pla1a.0.nr2)[1]
b <- range(pla1a.0.nr2)[2]
x <- tmp.sample1a[,3]
sex.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- tmp.sample1a[,4]
control.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- as.factor(tmp.sample1a[,5]) %>% as.numeric()
time.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type, 
               control.type,control.type, control.type, control.type, control.type,
               time.type, time.type, time.type, time.type, time.type)

image(
  cbind(pla1a.0.nr2[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]), heatmap(pla1a.0.nr2)$colInd ], 
        sidebar[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),]),
  col=redblue100,axes=F,main=glue("{TIS}2, f{P1}-HC, Sex-Control-Time"), asp =1)

plot(apply(pla1a.0.nr2[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),
     heatmap(pla1a.0.nr2)$colInd],1,median))

7.0.3 Transformations – only for the test samples. This can be revised to remove outlier samples

Strategy is not functionally programed (must revert log transformed back to linear gor pla1a.03a)

# pla1a.0.r3a<-pla1a.0.r3b<-pla1a.0.r3c<-pla1a.0.r3c2<-pla1a.0.r2b<-pla1a.0.r3d<-pla1a.0.r3d2<-pla1a.0.r2b2<-pla1a.0.r3d3 <-pla1a.0.2
pla1a.03a<-pla1a.03b<-pla1a.03c<-pla1a.03c2<-pla1a.02b<-pla1a.03d<-pla1a.03d2<-pla1a.02b2<-pla1a.03d3 <-pla1a.0.nr2
#pla1a.03c3<-pla1a.03b2<-pla1a.0.nr2

7.0.4 Examine the median and mean distributions (2)

tmp.s.median <- apply(pla1a.0.nr2,1, median)
tmp.s.mean <- apply(pla1a.0.nr2,1, mean)
plot(tmp.s.median,tmp.s.mean, asp = 1); abline(0,1)

7.0.5 Examine the median and sd distribution (1)

tmp.f.median <- apply(2^pla1a.0.nr2,2, median)
tmp.f.sd <- apply(2^pla1a.0.nr2,2, sd)
plot(y = tmp.f.sd, x = tmp.f.median,log="xy")

7.0.6 Examine the median and sd distribution (1; w/wo outlier samples)

tmp <- 2^pla1a.0.nr2[n.s,]
#dim(tmp)
tmp2.f.median <- apply(tmp,2, median)
tmp2.f.sd <- apply(tmp,2, sd)
plot(tmp.f.median,tmp2.f.median,log="xy")

plot(tmp.f.sd,tmp2.f.sd,log="xy")

7.1 Center by median, scale by standard deviation across features (3a)

for (i in 1:length(tmp.f.median)) {
  pla1a.03a[,i]<-(2^pla1a.0.nr2[,i]- tmp.f.median[i])/ tmp.f.sd[i]
}
plot(2^pla1a.0.nr2[,1],pla1a.03a[,1]) #spot-check, verified

7.1.1 N-P Heatmap: median-centered, sd-scaled imputed, run order (3a)

Notice the drift downward

# Run Order; Original HC
###########################
hmo <- heatmap(pla1a.0.nr2)$colInd

image(pla1a.03a[order(tmp.sample1a[,2]), hmo],
  col=redblue100,axes=F,main=glue("{TIS}3a, f{P1}-HC, Run-Order"), asp = 1)

# Run Order; New HC
###########################
hmo <- heatmap(pla1a.03a)$colInd

image(pla1a.03a[order(tmp.sample1a[,2]), hmo],
  col=redblue100,axes=F,main=glue("{TIS}3a, f{P1}-HC, Run-Order"), asp = 1)

plot(apply(pla1a.03a[order(tmp.sample1a[,2]),
     hmo],1,median)); abline(h = 0, col = 'blue')

7.1.2 Examine the median and sd distribution (2)

tmp.f.median <- apply(pla1a.0.nr2,2, median)
tmp.f.sd <- apply(pla1a.0.nr2,2, sd)
plot(tmp.f.median,tmp.f.sd, log = "xy")

7.1.3 Examine the median and sd distribution (2; w/wo outlier samples)

tmp <- pla1a.0.nr2[n.s,]
#dim(tmp)
tmp2.f.median <- apply(tmp,2, median)
tmp2.f.sd <- apply(tmp,2, sd)
plot(tmp.f.median,tmp2.f.median,log="xy")

plot(tmp.f.sd,tmp2.f.sd,log="xy")

7.2 Log2, Center by median, scale by standard deviation across features (3b)

for (i in 1:length(tmp.f.median)) {
  pla1a.03b[,i]<-(pla1a.0.nr2[,i]- tmp.f.median[i])/ tmp.f.sd[i]
}
plot(pla1a.0.nr2[,1],pla1a.03b[,1]) #verified

7.2.1 N-P Heatmap: median-centered, sd-scaled imputed, run order (3b)

# Run Order; Original HC
###########################
hmo <- heatmap(pla1a.0.nr2)$colInd

image(pla1a.03b[order(tmp.sample1a[,2]), hmo],
  col=redblue100,axes=F,main=glue("{TIS}3b, f{P1}-HC, Run-Order"), asp = 1)

# Run Order; New HC
###########################
hmo <- heatmap(pla1a.03b)$colInd

image(pla1a.03b[order(tmp.sample1a[,2]), hmo],
  col=redblue100,axes=F,main=glue("{TIS}3b, f{P1}-HC-NEW, Run-Order"), asp = 1)

plot(apply(pla1a.03b[order(tmp.sample1a[,2]),
     hmo],1,median))

7.2.2 N-P Heatmap: median-centered, sd-scaled imputed, Sex-Control-Time order (3b)

# Sex Control Time
a <- range(pla1a.03b)[1]
b <- range(pla1a.03b)[2]
x <- tmp.sample1a[,3]
sex.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- tmp.sample1a[,4]
control.type <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
x <- as.factor(tmp.sample1a[,5]) %>% as.numeric()
time.type  <- ( (b - a) * ((x - min(x))/((max(x) - min(x)))) + a) 
sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type, 
               control.type,control.type, control.type, control.type, control.type,
               time.type, time.type, time.type, time.type, time.type)

image(
  cbind(pla1a.03b[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]), heatmap(pla1a.0.nr2)$colInd ], 
        sidebar[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),]),
  col=redblue100,axes=F,main=glue("{TIS}3b, f{P1}-HC, Sex-Control-Time"), asp = 1)

plot(pla1a.03b[order(tmp.sample1a[,3],tmp.sample1a[,4],tmp.sample1a[,5]),])

7.2.3 Examine the sample median and sd distributions (3b)

Outlier sample has increased sd

boxplot(data.frame(t(pla1a.03b)))

plot(apply(pla1a.0.nr2,1,median), apply(pla1a.03b,1,median))

# remove outliers and determine if there is significant variablility in sample medians
glue("Outlier samples removed:")
## Outlier samples removed:
boxplot(data.frame(t(pla1a.03b[n.s,])))

7.2.4 TODO (in additional script): Collect sample medians independent of run-order outliers to determine if samples should be centered and/or scaled

7.2.5 TODO: Incorporate outlier sample removal

7.2.6 TODO: Incorporate flagged samples

7.2.7 TODO: Incorporate outlier feature removal from transformations

8 Save the processed data

8.1 Save the Data and Processing Decisions

# # log 2
# pla_rppos_pass1a_nr1 <- pla1a.0.nr1
# pla_rppos_pass1a_1 <- pla1a.0.1
# 
# # Test only
# #############
# # log2, imputed (test only)
# pla_rppos_pass1a_nr2 <- pla1a.0.nr2
# # linear, median-centered, sd-scaled
# pla_rppos_pass1a_3a <- pla1a.03a
# # log2, median-centered, sd-scaled
# pla_rppos_pass1a_3b <- pla1a.03b
feature_impute = names(feature_impute)

# The processing decisions
pla1a.0_vars <- data.frame(ds, site, tech, tis, NR, N, P, neg_vals, zero_vals, feature_na_filter,
P1, NR1, N1, paste(feature_impute, collapse = '; '), na_vals_impute, knn_k,
run_var, outlier_sample_n = length(o.s), outlier_samples =  paste(o.s, collapse = '; '), outlier_filter = o.f, 
internal_standards = paste(is, collapse = '; '),
comments = "No obvious sample outliers, data looks clean")

# visualize the processing decisions
if(knit_time){
  pla1a.0_vars %>%
    knitr::kable(format = "html") %>%
    scroll_box(width = "100%", height = "100%")
}
ds site tech tis NR N P neg_vals zero_vals feature_na_filter P1 NR1 N1 paste.feature_impute..collapse…….. na_vals_impute knn_k run_var outlier_sample_n outlier_samples outlier_filter internal_standards comments
pass1a umich rppos pla 101 77 170 0 0 10 170 101 77 0 10 0 1 90122013104 0.96 valine-13C5 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; methionine-13C5 [iSTD]; Car-D9 [iSTD]; phenylalanine-13C9 [iSTD]; arginine-13C6 [iSTD]; glucose-13C6 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Car(2:0)-D3 [iSTD]; Car(3:0)-D3 [iSTD]; Car(4:0)-D3 [iSTD]; Car(8:0)-D3 [iSTD]; fructose 1,6-bisphosphate-13C6 [iSTD]; Gibberellic acid [iSTD]; Car(16:0)-D3 [iSTD]; 24-epi-Brassinolide [iSTD] No obvious sample outliers, data looks clean 
# save the data
save(pla1a.0.nr1, pla1a.0.1, pla1a.0.nr2, pla1a.03a, pla1a.03b, pla1a.0_vars,
  file = glue("{WD}/data/UM-rppos/UM_rppos_processed_pla1a.0.0.RData"))

9 Session Info

warnings()
session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.2 (2020-06-22)
##  os       macOS  10.16                
##  system   x86_64, darwin17.0          
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Toronto             
##  date     2021-12-11                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source                                
##  assertthat     0.2.1    2019-03-21 [2] CRAN (R 4.0.2)                        
##  backports      1.2.1    2020-12-09 [2] CRAN (R 4.0.2)                        
##  broom          0.7.8    2021-06-24 [1] CRAN (R 4.0.2)                        
##  cachem         1.0.3    2021-02-04 [2] CRAN (R 4.0.2)                        
##  callr          3.7.0    2021-04-20 [1] CRAN (R 4.0.2)                        
##  cellranger     1.1.0    2016-07-27 [2] CRAN (R 4.0.2)                        
##  cli            3.0.1    2021-07-17 [1] CRAN (R 4.0.2)                        
##  coda           0.19-4   2020-09-30 [1] CRAN (R 4.0.2)                        
##  codetools      0.2-18   2020-11-04 [2] CRAN (R 4.0.2)                        
##  colorspace     2.0-2    2021-06-24 [1] CRAN (R 4.0.2)                        
##  crayon         1.4.1    2021-02-08 [2] CRAN (R 4.0.2)                        
##  curl           4.3.2    2021-06-23 [1] CRAN (R 4.0.2)                        
##  data.table     1.14.0   2021-02-21 [1] CRAN (R 4.0.2)                        
##  DBI            1.1.1    2021-01-15 [2] CRAN (R 4.0.2)                        
##  dbplyr         2.1.1    2021-04-06 [1] CRAN (R 4.0.2)                        
##  desc           1.3.0    2021-03-05 [1] CRAN (R 4.0.2)                        
##  devtools     * 2.4.2    2021-06-07 [1] CRAN (R 4.0.2)                        
##  digest         0.6.27   2020-10-24 [2] CRAN (R 4.0.2)                        
##  dplyr        * 1.0.7    2021-06-18 [1] CRAN (R 4.0.2)                        
##  ellipsis       0.3.2    2021-04-29 [1] CRAN (R 4.0.2)                        
##  evaluate       0.14     2019-05-28 [2] CRAN (R 4.0.1)                        
##  fansi          0.5.0    2021-05-25 [1] CRAN (R 4.0.2)                        
##  farver         2.1.0    2021-02-28 [1] CRAN (R 4.0.2)                        
##  fastmap        1.1.0    2021-01-25 [2] CRAN (R 4.0.2)                        
##  forcats      * 0.5.1    2021-01-27 [2] CRAN (R 4.0.2)                        
##  fs             1.5.0    2020-07-31 [2] CRAN (R 4.0.2)                        
##  generics       0.1.0    2020-10-31 [2] CRAN (R 4.0.2)                        
##  ggfittext      0.9.1    2021-01-30 [2] CRAN (R 4.0.2)                        
##  ggplot2      * 3.3.5    2021-06-25 [1] CRAN (R 4.0.2)                        
##  glue         * 1.4.2    2020-08-27 [2] CRAN (R 4.0.2)                        
##  gridExtra      2.3      2017-09-09 [2] CRAN (R 4.0.2)                        
##  gtable         0.3.0    2019-03-25 [2] CRAN (R 4.0.2)                        
##  haven          2.3.1    2020-06-01 [2] CRAN (R 4.0.2)                        
##  highr          0.9      2021-04-16 [1] CRAN (R 4.0.2)                        
##  hms            1.1.0    2021-05-17 [1] CRAN (R 4.0.2)                        
##  htmltools      0.5.1.1  2021-01-22 [2] CRAN (R 4.0.2)                        
##  httr           1.4.2    2020-07-20 [2] CRAN (R 4.0.2)                        
##  impute       * 1.62.0   2020-04-27 [1] Bioconductor                          
##  inline         0.3.19   2021-05-31 [1] CRAN (R 4.0.2)                        
##  inspectdf      0.0.11   2021-04-02 [1] CRAN (R 4.0.2)                        
##  jsonlite       1.7.2    2020-12-09 [2] CRAN (R 4.0.2)                        
##  kableExtra   * 1.3.1    2020-10-22 [2] CRAN (R 4.0.2)                        
##  knitr          1.33     2021-04-24 [1] CRAN (R 4.0.2)                        
##  labeling       0.4.2    2020-10-20 [2] CRAN (R 4.0.2)                        
##  lattice        0.20-41  2020-04-02 [2] CRAN (R 4.0.2)                        
##  lifecycle      1.0.0    2021-02-15 [1] CRAN (R 4.0.2)                        
##  loo            2.4.1    2020-12-09 [1] CRAN (R 4.0.2)                        
##  lubridate      1.7.10   2021-02-26 [1] CRAN (R 4.0.2)                        
##  magrittr       2.0.1    2020-11-17 [2] CRAN (R 4.0.2)                        
##  MASS           7.3-53   2020-09-09 [2] CRAN (R 4.0.2)                        
##  matrixStats    0.60.0   2021-07-26 [1] CRAN (R 4.0.2)                        
##  memoise        2.0.0    2021-01-26 [2] CRAN (R 4.0.2)                        
##  modelr         0.1.8    2020-05-19 [2] CRAN (R 4.0.2)                        
##  MotrpacBicQC * 0.5.2    2021-06-25 [1] Github (MoTrPAC/MotrpacBicQC@43fb293) 
##  munsell        0.5.0    2018-06-12 [2] CRAN (R 4.0.2)                        
##  mvtnorm        1.1-2    2021-06-07 [1] CRAN (R 4.0.2)                        
##  naniar         0.6.1    2021-05-14 [1] CRAN (R 4.0.2)                        
##  pillar         1.6.2    2021-07-29 [1] CRAN (R 4.0.2)                        
##  pkgbuild       1.2.0    2020-12-15 [2] CRAN (R 4.0.2)                        
##  pkgconfig      2.0.3    2019-09-22 [2] CRAN (R 4.0.2)                        
##  pkgload        1.2.1    2021-04-06 [1] CRAN (R 4.0.2)                        
##  plyr           1.8.6    2020-03-03 [2] CRAN (R 4.0.2)                        
##  prettyunits    1.1.1    2020-01-24 [2] CRAN (R 4.0.2)                        
##  processx       3.5.2    2021-04-30 [1] CRAN (R 4.0.2)                        
##  progress       1.2.2    2019-05-16 [2] CRAN (R 4.0.2)                        
##  ps             1.6.0    2021-02-28 [1] CRAN (R 4.0.2)                        
##  purrr        * 0.3.4    2020-04-17 [2] CRAN (R 4.0.2)                        
##  R6             2.5.0    2020-10-28 [2] CRAN (R 4.0.2)                        
##  Rcpp           1.0.7    2021-07-07 [1] CRAN (R 4.0.2)                        
##  RcppParallel   5.1.4    2021-05-04 [1] CRAN (R 4.0.2)                        
##  readr        * 1.4.0    2020-10-05 [2] CRAN (R 4.0.2)                        
##  readxl         1.3.1    2019-03-13 [2] CRAN (R 4.0.2)                        
##  remotes        2.4.0    2021-06-02 [1] CRAN (R 4.0.2)                        
##  reprex         2.0.0    2021-04-02 [1] CRAN (R 4.0.2)                        
##  reshape2       1.4.4    2020-04-09 [2] CRAN (R 4.0.2)                        
##  rethinking   * 2.13     2021-08-13 [1] Github (rmcelreath/rethinking@2acf2fd)
##  rlang          0.4.11   2021-04-30 [1] CRAN (R 4.0.2)                        
##  rmarkdown      2.6      2020-12-14 [2] CRAN (R 4.0.2)                        
##  rprojroot      2.0.2    2020-11-15 [2] CRAN (R 4.0.2)                        
##  rstan        * 2.21.2   2020-07-27 [1] CRAN (R 4.0.2)                        
##  rstudioapi     0.13     2020-11-12 [2] CRAN (R 4.0.2)                        
##  rvest          1.0.0    2021-03-09 [1] CRAN (R 4.0.2)                        
##  scales         1.1.1    2020-05-11 [2] CRAN (R 4.0.2)                        
##  sessioninfo    1.1.1    2018-11-05 [2] CRAN (R 4.0.2)                        
##  shape          1.4.6    2021-05-19 [1] CRAN (R 4.0.2)                        
##  StanHeaders  * 2.21.0-7 2020-12-17 [1] CRAN (R 4.0.2)                        
##  stringi        1.7.3    2021-07-16 [1] CRAN (R 4.0.2)                        
##  stringr      * 1.4.0    2019-02-10 [2] CRAN (R 4.0.2)                        
##  testthat       3.0.4    2021-07-01 [1] CRAN (R 4.0.2)                        
##  tibble       * 3.1.3    2021-07-23 [1] CRAN (R 4.0.2)                        
##  tidyr        * 1.1.3    2021-03-03 [1] CRAN (R 4.0.2)                        
##  tidyselect     1.1.1    2021-04-30 [1] CRAN (R 4.0.2)                        
##  tidyverse    * 1.3.1    2021-04-15 [1] CRAN (R 4.0.2)                        
##  usethis      * 2.0.1    2021-02-10 [1] CRAN (R 4.0.2)                        
##  utf8           1.2.2    2021-07-24 [1] CRAN (R 4.0.2)                        
##  V8             3.4.2    2021-05-01 [1] CRAN (R 4.0.2)                        
##  vctrs          0.3.8    2021-04-29 [1] CRAN (R 4.0.2)                        
##  viridisLite    0.4.0    2021-04-13 [1] CRAN (R 4.0.2)                        
##  visdat         0.5.3    2019-02-15 [2] CRAN (R 4.0.2)                        
##  webshot        0.5.2    2019-11-22 [2] CRAN (R 4.0.2)                        
##  withr          2.4.2    2021-04-18 [1] CRAN (R 4.0.2)                        
##  xfun           0.24     2021-06-15 [1] CRAN (R 4.0.2)                        
##  xml2           1.3.2    2020-04-23 [2] CRAN (R 4.0.2)                        
##  yaml           2.2.1    2020-02-01 [2] CRAN (R 4.0.2)                        
## 
## [1] /Users/Alec/Library/R/4.0/library
## [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library